data {
  m.1 <- m - 1
}
model {
  for (i in 1:n) {
    y[i, 1:m] ~ dmulti(theta[i, 1:m], k[i])
    for (j in 1:m.1) {
      theta[i, j] <- exp.theta.star[i, j] / sum.exp.theta.star[i]
    }
    theta[i, m] <- 1 / sum.exp.theta.star[i]
    sum.exp.theta.star[i] <- sum(exp.theta.star[i, 1:m.1]) + 1
    for (j in 1:m.1) {
      exp.theta.star[i, j] <- exp(theta.star[i, j])
    }
  }
  for (i in 2:n) {
    # theta.star[i, 1:m.1] ~ dmnorm(mu[i, 1:m.1], Omega.star)
    for (j in 1:m.1) {
      theta.star[i, j] ~ dnorm(mu[i, j], Omega.star[j])
      # mu[i, j] <- psi[j] + phi[j] * theta.star[i - 1, j] + inprod(beta[j, ], x[i, ])
      mu[i, j] <- theta.star[i - 1, j] + inprod(beta[j, ], x[i, ])
    }
  }
  theta.star[1, 1:m.1] ~ dmnorm(mu0, Omega.star0)
  for (i in 1:m.1) {
    for (j in 1:l) {
      beta[i, j] ~ dnorm(0, 0.01)
    }
    # psi[i] ~ dnorm(0, 0.0001)
    # phi[i] ~ dnorm(0, 0.0001)
    Omega.star[i] ~ dscaled.gamma(100, 2)
    Omega[i] <- pow(Omega.star[i], -2)
  }
  # Omega.star ~ dscaled.wishart(rep(100, m.1), 2)
  # Omega <- inverse(Omega.star)
}